clc;
clear;

h = 0.5;
left = -1;
right = 1;
bottom = 0;
top = 1;
Nx = (right-left)/h;
Ny = (top-bottom)/h;

basis_type = "P2";

%% mesh_read
filename = 'mesh_hole.msh';
fileID = fopen(filename, 'r');
data = textscan(fileID, '%f %f %f %f', 'Delimiter', ' ', 'MultipleDelimsAsOne', true);
fclose(fileID);
Np = data{1}(1);
Ne = data{2}(1);
node1 = [data{1}(2:Np+1), data{2}(2:Np+1)];
elem1 = [data{1}(Np+2:Np+Ne+1), data{2}(Np+2:Np+Ne+1), data{3}(Np+2:Np+Ne+1)];
bdary = [data{1}(Np+Ne+2:end), data{2}(Np+Ne+2:end)];
bde = boundaryedges(elem1, bdary);
[node2, elem2, bdn] = process_grid_interpolation(node1, elem1, bde, basis_type);
Ne = size(elem2,1);
Nb = size(node2,1);
Nlb = size(elem2,2);

% elem2 =
%      8    10     9    11
%      3     7     4    12
%      4     8     5    13
%      4     2     1    14
%      5     2     4    15
%      4     1     3    16
%     10     8     7    17
%      8     6     5    18
%      7     8     4    19
%      9     6     8    20

% elem2 =
%      8    10     9    11    12    13
%      3     7     4    14    15    16
%      4     8     5    17    18    19
%      4     2     1    20    21    22
%      5     2     4    23    20    19
%      4     1     3    22    24    16
%     10     8     7    11    25    26
%      8     6     5    27    28    18
%      7     8     4    25    17    15
%      9     6     8    29    27    13

% node1 =
%    -0.0000   -1.0000
%    -0.7071   -0.7071
%     0.7071   -0.7071
%    -0.0000   -0.3333
%    -1.0000    0.0000
%    -0.7071    0.7071
%     1.0000         0
%    -0.0000    0.3333
%     0.0000    1.0000
%     0.7071    0.7071

% node2 =
%    -0.0000   -1.0000
%    -0.7071   -0.7071
%     0.7071   -0.7071
%    -0.0000   -0.3333
%    -1.0000    0.0000
%    -0.7071    0.7071
%     1.0000         0
%    -0.0000    0.3333
%     0.0000    1.0000
%     0.7071    0.7071
%     0.3536    0.5202
%     0.3536    0.8536
%    -0.0000    0.6667
%     0.8536   -0.3536
%     0.5000   -0.1667
%     0.3536   -0.5202
%    -0.0000   -0.0000
%    -0.5000    0.1667
%    -0.5000   -0.1667
%    -0.3536   -0.5202
%    -0.3536   -0.8536
%    -0.0000   -0.6667
%    -0.8536   -0.3536
%     0.3536   -0.8536
%     0.5000    0.1667
%     0.8536    0.3536
%    -0.3536    0.5202
%    -0.8536    0.3536
%    -0.3536    0.8536

if basis_type == "P1b"
    mod = [1, 2, 4;
        2, 3, 4;
        3, 1, 4];
    elem_p1b = reshape(elem2(:, mod(:)), 3*Ne, 3);
elseif basis_type == "P2"
    mod = [1, 4, 6;
        2, 5, 4;
        3, 6, 5;
        4, 5, 6];
    elem_p2 = reshape(elem2(:, mod(:)), 4*Ne, 3);
end

pde = process_pde(1);
mesh.basis_type = basis_type;
mesh.node = node2;
mesh.elem = elem2;
mesh.Ne = Ne;
mesh.Np = Np;
mesh.Nb = Nb;
process_showmesh(mesh);

figure;
x = node2(:,1);
y = node2(:,2);
U = pde.u1(x, y);
V = pde.u2(x, y);
P = pde.p(node1(:,1), node1(:,2));
if basis_type == "P1"
    trisurf(elem2, x, y, U);
    figure;
    trisurf(elem2, x, y, V);
elseif basis_type == "P1b"
    trisurf(elem_p1b, x, y, U);
    figure;
    trisurf(elem_p1b, x, y, V);
elseif basis_type == "P2"
    trisurf(elem_p2, x, y, U);
    figure;
    trisurf(elem_p2, x, y, V);
end
figure;
trisurf(elem1, node1(:,1), node1(:,2), P);



% [X, Y] = meshgrid(linspace(min(x), max(x), size(x,1)), linspace(min(y), max(y), size(y,1)));
% in_poly = inpolygon(X(:), Y(:), node1(bdary(:),1), node1(bdary(:),2));
% in_poly = reshape(in_poly, size(X));
% F_Ux = scatteredInterpolant(x, y, U);
% F_Uy = scatteredInterpolant(x, y, V);
% F_P = scatteredInterpolant(node1(:,1), node1(:,2), P);
% Ux = F_Ux(X, Y);
% Uy = F_Uy(X, Y);
% Ped = F_P(X, Y);
% Ux(~in_poly) = NaN;
% Uy(~in_poly) = NaN;
% Ped(~in_poly) = NaN;


% figure;
% streamslice(X, Y, Ux, Uy);
% title('流线图');

% figure;
% quiver(X, Y, Ux, Uy);
% title('速度矢量场');

% figure;
% contourf(X, Y, Ped);
% title('压力场等值线图');

% figure;
% surf(X, Y, Ux);

% figure;
% contourf(X, Y, Ux);






%% process_pde
function PDE = process_pde(index)
switch index
    case 1
        a = 1;
        PDE.index = index;
        PDE.nu = a;
        PDE.u1 = @(x,y) x.^2 .* y.^2 + exp(-y);
        PDE.u1dx = @(x,y) 2*x.*y.^2;
        PDE.u1dy = @(x,y) 2*y.*x.^2 - exp(-y);
        PDE.u2 = @(x,y) -2/3*x.*y.^3 + 2 - pi*sin(pi*x);
        PDE.u2dx = @(x,y) - pi^2*cos(pi*x) - (2*y.^3)/3;
        PDE.u2dy = @(x,y) -2*x.*y.^2;
        PDE.p = @(x,y) -(2-pi*sin(pi*x)).*cos(2*pi*y);
        PDE.pdx = @(x,y) pi^2*cos(pi*x).*cos(2*pi*y);
        PDE.pdy = @(x,y) -2*pi*sin(2*pi*y).*(pi*sin(pi*x) - 2);
        PDE.f1 = @(x,y) -2*a*x.^2 - 2*a*y.^2 - a*exp(-y) + pi^2*cos(pi*x).*cos(2*pi*y);
        PDE.f2 = @(x,y) 4*a*x.*y - a*pi^3*sin(pi*x) + 2*pi*(2-pi*sin(pi*x)).*sin(2*pi*y);
        PDE.domain = [0, 1, -0.25, 0];
        PDE.bdary = [0, 0, 0, 0];
    case 2
        a = 1;
        PDE.index = index;
        PDE.nu = a;
        PDE.u1 = @(x,y) x.^2 .* y.^2 + exp(-y);
        PDE.u1dx = @(x,y) 2*x.*y.^2;
        PDE.u1dy = @(x,y) 2*y.*x.^2 - exp(-y);
        PDE.u2 = @(x,y) -2/3*x.*y.^3 + 2 - pi*sin(pi*x);
        PDE.u2dx = @(x,y) - pi^2*cos(pi*x) - (2*y.^3)/3;
        PDE.u2dy = @(x,y) -2*x.*y.^2;
        PDE.p = @(x,y) -(2-pi*sin(pi*x)).*cos(2*pi*y);
        PDE.pdx = @(x,y) pi^2*cos(pi*x).*cos(2*pi*y);
        PDE.pdy = @(x,y) -2*pi*sin(2*pi*y).*(pi*sin(pi*x) - 2);
        PDE.f1 = @(x,y) -2*a*x.^2 - 2*a*y.^2 - a*exp(-y) + pi^2*cos(pi*x).*cos(2*pi*y);
        PDE.f2 = @(x,y) 4*a*x.*y - a*pi^3*sin(pi*x) + 2*pi*(2-pi*sin(pi*x)).*sin(2*pi*y);
        PDE.domain = [0, 1, -0.25, 0];
        PDE.bdary = [0, 1, 0, 0];
    otherwise
        error("Invalid PDE index.");
end
fprintf("PDE index: %d\n", PDE.index);
fprintf("PDE domain: xmin,xmax,ymin,ymax\n");
disp(PDE.domain);
fprintf("PDE boundary: bottom,right,top,left (0=Dirichlet,1=Neumann,2=Robin)\n");
disp(PDE.bdary);
end

%% process_showmesh
function process_showmesh(mesh)
node = mesh.node;
elem = mesh.elem;
figure;
% mesh
h = patch('Faces', elem(:,1:3), 'Vertices', node(1:mesh.Nb,:));
set(h, 'facecolor', 'w', 'edgecolor', 'k');
axis equal;
axis tight;
hold on;
% elem
center = zeros(size(elem,1),2);
for i = 1:size(elem,1)
    V = node(elem(i,:),:);
    center(i,:) = polygon_centroid(V);
end
text(center(:,1), center(:,2), int2str((1:mesh.Ne)'), 'FontSize', 10, 'FontWeight', 'bold', 'Color', 'r');
% node
shift = 0.001;
range = (1:mesh.Np)';
plot(node(range,1), node(range,2), 'k.', 'MarkerSize', 10);
text(node(range,1)+shift, node(range,2)+shift, int2str(range), 'FontSize', 10, 'FontWeight', 'bold', 'Color', 'k');
if mesh.basis_type == "P1b" || mesh.basis_type == "P2"
    plot(node(mesh.Np+1:end,1), node(mesh.Np+1:end,2), 'b.', 'MarkerSize', 10);
    text(node(mesh.Np+1:end,1)+shift, node(mesh.Np+1:end,2)+shift, int2str((mesh.Np+1:mesh.Nb)'), 'FontSize', 10, 'FontWeight', 'bold', 'Color', 'b');
end
hold off
end








%% process_square
function [node, elem, bde] = process_square(left, right, bottom, top, Nx, Ny)
[X, Y] = meshgrid(linspace(left, right, Nx+1), linspace(bottom, top, Ny+1));
node = [X(:), Y(:)];
Np = size(node,1);
k = (1:Np-Ny-1)';
cut = (Ny+1)*(1:Nx);
k(cut) = [];
elem = [k, k+Ny+1, k+1;
    k+1, k+Ny+1, k+Ny+2];
elem = sortrows(elem);
Ne = size(elem,1);
bd_bottom = [1:2*Ny:Ne;
    1:Ny+1:Np-Ny-1;
    Ny+2:Ny+1:Np]';
bd_right = [(Nx-1)*2*Ny+2:2:Ne;
    Np-Ny:Np-1;
    Np-Ny+1:Np]';
bd_top = [2*Ny:2*Ny:Ne;
    2*Ny+2:Ny+1:Np;
    Ny+1:Ny+1:Np-1]';
bd_top = sortrows(bd_top,1,"descend");
bd_left = [1:2:2*Ny-1;
    2:Ny+1;
    1:Ny]';
bd_left = sortrows(bd_left,1,"descend");
bde = [bd_bottom; bd_right; bd_top; bd_left];
end

%% process_grid_interpolation
function [new_node, new_elem, new_bdn] = process_grid_interpolation(node, elem, bde, basis_type)
bdn = unique(bde(:,2:3), "stable");
switch basis_type
    case "P1"
        new_node = node;
        new_elem = elem;
        new_bdn = bdn;
    case "P1b"
        P = zeros(size(elem,1),2);
        for k = 1:size(elem,1)
            V = [node(elem(k,:),:)];
            P(k,:) = polygon_centroid(V);
        end
        new_col = (size(node,1)+1:size(node,1)+size(elem,1))';
        new_node = [node; P];
        new_elem = [elem, new_col];
        new_bdn = bdn;
    case "P2"
        v0 = reshape(elem', 3*size(elem,1), 1);
        v1 = reshape(([elem(:,2:3), elem(:,1)])', 3*size(elem,1), 1);
        totalEdge = [v0,v1];
        index1 = zeros(size(totalEdge,1),2);
        k = 1;
        for i = 1:size(totalEdge,1)
            if ~all(totalEdge(i,:) == [0, 0])
                for j = i+1:size(totalEdge,1)
                    if all(totalEdge(i,:) == flip(totalEdge(j,:)))
                        totalEdge(j,:) = [0, 0];
                        index1(k,:) = [j, i];
                        k = k+1;
                    end
                end
            end
        end
        uniqueEdge = totalEdge(~(totalEdge(:,1) == 0),:);
        P = (node(uniqueEdge(:,1),:) + node(uniqueEdge(:,2),:))/2;
        index0 = setdiff((1:size(totalEdge,1))', sort(index1(:)));
        index1 = sortrows(index1(~(index1(:,1) == 0),:));
        index2 = [find(totalEdge(:,1) ~= 0), (1:size(P,1))'];
        new_col = zeros(size(elem,1), 3);
        k = 1;
        for i = 1:3:size(totalEdge,1)-2
            new_col(k,:) = [i, i+1, i+2];
            k = k+1;
        end
        for k = 1:size(index1,1)
            idx = (new_col == index1(k,1));
            new_col(idx) = index1(k,2);
        end
        for k = 1:size(index2,1)
            idx = (new_col == index2(k,1));
            new_col(idx) = index2(k,2);
            idx = (index0 == index2(k,1));
            index0(idx) = index2(k,2);
        end
        [~,index3] = sortrows(bde(:,1));
        [~,index3] = sortrows(index3);
        index0 = index0(index3);
        index0 = index0 + size(node,1);
        new_node = [node; P];
        new_elem = [elem, new_col+size(node,1)];
        new_bdn = reshape([bdn'; index0'],[],1);
    otherwise
        error("Invalid basis type.");
end
end

%% polygon_centroid
function center = polygon_centroid(V)
n = size(V,1);
A = 0;
center_x = 0;
center_y = 0;
for i = 1:n
    j = mod(i,n) + 1;
    xi = V(i,1);
    yi = V(i,2);
    xj = V(j,1);
    yj = V(j,2);
    factor = xi * yj - xj * yi;
    A = A + factor;
    center_x = center_x + (xi + xj) * factor;
    center_y = center_y + (yi + yj) * factor;
end
center_x = center_x / (3 * A);
center_y = center_y / (3 * A);
center = [center_x, center_y];
end

%% boundaryedge
function bde = boundaryedges(elem, bdary)
E_circle = [elem, elem(:,1)];
edges = [E_circle(:,1:2); E_circle(:,2:3); E_circle(:,3:4)];
[~, loc] = ismember(edges, bdary, 'rows');
E_bd = reshape(loc, size(elem,1), 3);
[rows, cols] = find(E_bd);
index0 = zeros(size(rows,1), 1);
for i = 1:size(rows,1)
    index0(E_bd(rows(i), cols(i))) = rows(i);
end
bde = [index0, bdary];
end
